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Inclusive cross sections for J/^p production in proton- proton collisions were calculated in the 
fct-factorization approach for the RHIC energy. Several mechanisms were considered, including 
direct color-singlet mechanism, radiative decays of Xc mesons, decays of tp' , open-charm associated 
production of J/ip as well as weak decays of B mesons. Different unintegrated gluon distributions 
from the literature were used. We find that radiative Xc decays and direct color-singlet contributions 
constitute the dominant mechanism of J/^ production. These process cannot be consistently treated 
within coUinear-factorization approach. The results are compared with recent RHIC data. The new 
precise data at small transverse momenta impose stringent constraints on UGDFs. Some UGDFs 
Ch , are inconsistent with the new data. The Kwieciriski UGDFs give the best description of the data. In 

, ^ ■ order to verify the mechanism suggested here we propose J/tp - jet correlation measurement and an 

independent measurement of Xc meson production in tt^tt" and/or K'^ K~ decay channels. Finally, 
we address the issue of J/ip spin alignment. 

■ PACS numbers: 12.38Bx,13.85Ni,14.40Gx 

^\ 

CD I. INTRODUCTION 

For the last decade, the inclusive production of J/tp mesons was a serious theoretical puzzle challenging our un- 
O^l ' derstanding of QCD, parton model, and the bound state formation dynamics. The roots of the puzzling J /ip history 
', trace back to the middle 1990s, when the data on J and T hadroproduction cross sections revealed a more 

■ than one order-of-magnitude discrepancy with theoretical expectations. This fact has induced extensive theoretical 
activity and led to the introduction of new production mechanisms, known as the color-octet model [1, Q and gluon 
vector dominance model Since then, the color-octet model has been believed to give the most likely explanation 
of the quarkonium production phenomena, although there were also some indications that it was not working well. 
The situation became even more intriguing after the measurements of J /ip spin alignment 01" H have been carried 
out showing inconsistency with the newly accepted theory. 

, At the same time, it has been shown that the incorporation of the usual color-singlet production scheme with 
1^ ' the fcj -factorization approach can provide a reasonable and consistent picture of the phenomenon under study in its 
.5^ ] entirety. Whithin the latter approach, a good description of data on the production of J/'0, Xcj and T mesons both 
. at the Tevatron p, jioj and HERA [ll| has been achieved, and even a solution to the J/ tp spin alignment problem has 
?H ■ been guessed P. |l2|. The issue of the quarkonium production mechanism continues to be under intense debate. 

Recently, the PHENIX collaboration at the BNL Relativistic Heavy Ion Collider (RHIC) has measured inclusive 
J/ip production in elementary proton-proton collisions [l3| ■ While for the RHIC community the elementary pp cross 
section is only the baseline for the nuclear case, we wish to demonstrate that the elementary data by itself constitute 
a very valueable information about QCD dynamics in the region of intermediate x ~ 10~^-10~^. In our paper we 
present a detailed analysis of RHIC data based on the fct-factorization approach and a large variety of unintegrated 
gluon distribution functions (UGDFs). We show that the new precise data at small J/ip transverse momenta impose 
stringent constraints on UGDFs and, consequently, stimulate better understanding of the underlying gluon dynamics. 

The outline of the paper is the following. In Sec. II we describe the production mechanisms employed in our 
analysis and discuss the different parametrizations of UGDFs. In Sec. HI we compare our theoretical predictions with 
experimental results and derive new predictions on the quantities which at yet have not been measured but could 
serve as important cross-check of our understanding of the reaction mechanism. Our findings and recommendations 
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for the forthcoming experiments are summarized in Sec. IV. 



II. FORMALISM 



A. Different mechanisms of J/t/j production 

In this paper, we take into account a number of different mechanisms leading to the appearance of J /if) mesons in 
the final state (of course, they are not thought to be all of equal importance). The considered mechanisms are the 
following. 

Direct color-singlet J/?/; production via gluon-giuon fusion 

g + g^J/i) + g] (1) 

direct production of ip' meson 

g + g^i^' + g (2) 

and its subsequent decay ip' J/ip + X; 

production of P-wave charmonium states Xcj (^^ 0, 1, 2) 



9 + 9^ XcJ (3) 



followed by their radiative decays Xcj —> J/'4' + 7; 
production of b quarks and antiquarks 



g + g^b + b (4) 

followed by their fragmentation into B mesons and subsequent weak decays B — > J/^p + X; 
production of J/ijj mesons in association with unbound charmed quarks 

g + g^ J/ip + C + C. (5) 

Examples of the relevant Feynman diagrams for all the mentioned processes are shown in Fig. 1. Every subprocess 
is accompanied by the emission of gluon jets, as is shown in Fig. 2. 

In general, there could also exist color-octet contributions. The latter cannot be calculated from the first principles 
and are usually estimated from fits to existing data. It has been shown already that, within the fej-factorization 
approach, these contribuions are consistent with zero both at the Tevatron and HERA In view of the 

uncertainties coming from other contributions we find it not useful to include color-octet contributions in the present 
analysis. 

A few words are in order to describe the formation of cc bound states. First of all, it should be noted that the 
amplitudes of the subprocesses dJ)-®, (O contain projection operators J{S, L), which guarantee the proper quantum 
numbers of the cc state under consideration. These operators read for the different spin and orbital angular momentum 
states [11 Hi: 

-/5i^)c + mc)/mlil'^ , (6) 
4iS.){^, + m,)/m'J\ (7) 
- ^c) 4iSz) {i>c + m,)/w?J^ , (8) 

where is the mass of the specifically considered cc state and TOc = the mass of the charmed quark (always 

set equal to 1/2 of the meson mass, as is required by the nonrelativistic bound-state model). 

States with various projections of the spin momentum onto the z axis are represented by the polarization vector 

The probability for the two quarks to form a meson depends on the bound state wave function '^{q). In the 
nonrelativistic approximation which we are using here, the relative momentum q of the quarks in the bound state is 
treated as a small quantity. So, it is useful to represent the quark momenta as Pc = Vi'l^ + 9i Vc— Pv/^ ~ 1- Then, 
we multiply the matrix elements by ^"((7) and perform integration with respect to q. The integration is performed 
after expanding the integrand around 5 = 0: 

M{q) = A^|,=o + (9X/ag")|,=o9" + ■ ■ ■ • (9) 



J('5o) 


EE J(5= 


=0,L= 


=0) 


J(351) 


= J(5= 


=1,L= 


=0) 




^ ,7(5= 


=1,L= 


=1) 





d) 





e) 



FIG. 1: Processes included in our approach: (a) direct color-singlet production, (b) production of Xc mesons, (c) open bottom 
quark production, (d) open-charm associated production, (e) color-octet production 



Since the expressions for Ai\q^07 {9.M/dq°')\q=Q, etc. are no longer dependent on q, they may be factored outside the 
integral sign. A term- by-term integration of this series then yields [l5| : 



_fq_ 

(27r)3 



7^(x = 0), 



V3 



= -ie"(L,) ^ n'{x = 0), 



An 



(10) 
(11) 



etc., where TZ{x) is the radial wave function in the coordinate representation (the Fourier transform of '^{q)). The 
first term contributes only to S"- waves, but vanishes for P- waves because TZp{0) = 0. On the contrary, the second 
term contributes only to P- waves, but vanishes for 5- waves because Ti'g{0) = 0. States with various projections of 
the orbital angular momentum onto the z axis are represented by the polarization vector e{Lz). The numerical values 
of the wave functions are either known from the leptonic decay widths (for J/ip and ip' mesons) or can be taken 
from potential models (for Xcj mesons). Including radiative corrections changes the values of the wave functions by 
a factor of 2 (the NLO result compared to the LO result), and so, one can also expect large effect from higher order 
corrections. This leads to a sizeable theoretical uncertainty, which, on the other hand, can only affect the absolute 
normalization but not the shape of the pt spectrum. 

When calculating the spin average of the matrix elements squared, we adopt the fcf-factorization prescription [l6j 



FIG. 2: Application of UGDFs to inclusive production of J/i/; (left) and Xc (right). The upper and the lower parts of these 
diagrams are included in the kt evolution of gluon densities. The emitted gluons can realise in the final state hadronic jets. 



for the off-shell gluon spin density matrix: 

W^^k'^KI\kt\'', (12) 

where kt is the component of the gluon momentum perpendicular to the beam axis, and the bar stands for the averaging 
over the gluon spin. In the coUinear limit, when fct — > 0, this expression converges to the ordinary e^e*'' = — ^ g^". 
In all other respects, the evaluation of the diagrams is straightforward and follows the standard QCD Feynman rules. 
This has been done using the algebraic manipulation systems FORM [l3] and REDUCE [Tsj . 
For the direct production mechanism ^ the fully differential cross section reads 



7 I /VI \ (HI — luii i I ^ 

64 

spina colors 

•KTg{xi,klt, y?) Tg{x2,klt, y?) dkltdkltdp^^dyzdy^'^'^'^ , (13) 

ZTT ZVr ZTT 

where 4>i and ^2 are the azimuthal angles of the initial gluons, and and 0^ the rapidity and the azimuthal angle of 
J/il) particle. The explicit expressions for the parton level matrix elements \M.{gg "^ff)!^ are presented in Ref. 
The phase space physical boundary is determined by the inequality [l^ 

G{s,iklkf,klml)<0, (14) 

with ki, k2 and fcs being the initial and final gluon momenta, s — {ki + fc2)^, t — {ki — p^)^, and G is the 
standard kinematic function [T^ . The initial gluon momentum fractions xi and X2 appearing in the unintegrated 
gluon distribution functions J^g{xi, kf t, fi^) are calculated from the energy-momentum conservation laws in the light 
cone projections: 

{ki + k2)E+pu = xi^/s = m^Texp(y^) -|- |A:3t| exp(y3), 

(15) 

{ki + k2)E-p^^ = X2VS = m^Texp(-?;^) -f IfcstI exp(-y3), 

where m^r = (m-^ + Ip^prl'^y^'^ ■ 

The production scheme of ip' meson ^ is identical to that of J/ip , and only the numerical value of the wave 
function |7?,(0)P is different. In both cases, the values of the wave functions were extracted from the known leptonic 
decay widths [23| using the formula |7^(0)p = Teem'^/iia'^el) [1 - 16Q;s/(37r)] and were set equal to |7^J/^(0)p = 0.8 
GeV'^ for J/tp meson, and |7?.0'(O)p ~ 0.4 GeV'^ for tp' meson. To calculate the feed-down to J/ip states, the ip' 
production cross section has to be multiplied by the branching fraction Br{%p' J/ijjX) — 56% [2(j |. 
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For the production of Xcj mesons via the subprocess ([3]) we have 

^^ 127r2a2|7^'(0)|2 i i ^ 

xTg{xi,klt, Tg{x2,klt, dkjtdkltdy^^^. (16) 

The squares of the matrix elements, as being too lengthy, are not presented there but the full fortran code is available 
on request. The numerical value of the wave function is taken from the potential model [2l|: |7^^(0)p = 0.075 GeV^ 
The decay branchings to J/tp meson are known to be f^oj Br{xcj J/i'l) = 0.013, 0.35, and 0.20 for J = 0, 1, and 2, 
respectively. Here, the off-shell gluon flux factor is defined as F = 2A^/^(s, kf, according to the general definition 
given by Eq.(2.3) in Ref. [l^. For all other subprocesses one can use the approximations A^/^(s, kf, fc|) ~ s ~ X1X2S, 
but they are not suitable for the present case because the invariant mass of the final state is small and the difference 
between s = and X1X2S = rn'^^ = + Pt can make pronounced effect on the px spectrum. The numerical 
accuracy of the above definition was tested in a toy calculation regarding the leptonic production of Xcj mesons via 
photon-photon fusion: e + e ^ e' + e' + Xc- We have compared the exact 0{a'^) result with a number of calculations 
based on Equivalent Photon Approximation and using different definitions of the effective photon flux (such as = 2s, 
F = 2xiX2S, etc.). We find that the "A^/^" definition is in the best agreement with exact calculation. 
For the production of beauty quarks in ^ wc have 

daipp ^ bbX) = 1^ i ^ 1 ^ |A1 (53 -> bb)f 

■xTg{xl,kl^, Tg{x2,kl^, dkl^ dkjt dpfj. dyt dy^ ^ ^ ^ . (17) 

The explicit expressions for the parton level matrix elements \A4{gg — > bb)\^ can be found elsewhere [l^]. In calcu- 
lations the 6-quark mass was set to rub — 4.5 GeV. Further on, the produced 6-quarks undergo fragmentation into 
B-mesons according to the Peterson fragmentation function [23j with e=0.006. The outgoing B-mesons undergo then 
a decay according to the three body decay mode B ^ J/^ + K + tt, to which the net effective branching fraction 
[20t was attributed: Br{b J/i'X) = 1.16% (resp., Br{b X) = 0.48%). This decay mode was taken as a 

typical representative for all _B-meson decays. As the decay matrix elements are unknown, the decays were generated 
according to the phase space. However, the fine details of fragmentation and decay are rather unimportant for our 
purposes, because 5-quarks play only marginal role at RHIC energies, except large transverse momenta of J /ip or ip' . 
We shall discuss the region of the large transverse momenta somewhat later. 
Finally, for the charm-associated production ([5]) wc write 

da{pp ^ccX) ^ ^ |7^(0)pi E ^ E 1-^(5.9 - V'cc)^ 

spins colors 

2 2^ J 2 2^ 17 2 77 2 J 2 7 2 7 ^"^1 d<p2 dtf)^ d(j)c 

xTg{xi,k^t,^j. )Tg[x2,k2t,tJ. ) dkudk2tdp^T<iPcTdy^dycdys- — t;—-^- (18) 

ZTT ZTT ZTT ZTT 

The explicit expressions for the parton level matrix elements \Ai{gg ipcc)\'^ as well as detailed description of the 
kinematics are presented in Ref. [23 |. 

To close the description of the production mechanisms, we wish to state that we do not consider explicitly color- 
octet contributions in the present analysis. In fact, we know no data which would clearly manifest the presence of 
color-octet contributions. On the contrary, the numerical fits of the color-octet matrix elements based on the Tevatron 
and HERA data are incompatible with each other. Moreover, a conflict between the model predictions and the data 
on J/ip spin alignment indicate that the production of vector quarkonia is certainly not dominated by the color-octet 
mechanism. Some small contribution is not excluded but cannot be calculated from first principles. 



B. Unintegrated gluon distributions 

In general, there are no simple relations between unintegrated and integrated parton distributions. Some of UGDFs 
in the literature are obtained based on familiar coUinear distributions, some are obtained by solving evolution equa- 
tions, some are just modelled or some are even parametrized. A brief review of unintegrated gluon distributions 
(UGDFs) that will be used also here can be found in Ref. ^] . 
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At very low x the unintegrated gluon distributions are believed to fulfil BFKL equation '26'|. Here in the practical 
applications we shall use a simple parametrization \27\ for the numerical solution i28 j | and use acronym BFKL. Another 
distribution closely related to the BFKL approach was constructed by Bliimlein [23| . 

At large energies (small x) one expects in addition saturation effects due to gluon recombinations. A simple 
parametrization of unintegrated gluon distribution in the proton can be obtained based on the Golec-Biernat-Wiisthoff 
parametrization of the dipole-nucleon cross section with parameters fitted to the HERA data. The dipole-nucleon 
cross section can be transformed to corresponding unintegrated gluon distribution. The resulting gluon distribution 
can be found in [s^l- In the following we call it GBW UGDF for brevity. Another parametrization, also based on 
the idea of gluon saturation, was proposed in [3l|. In contrast to GBW approach [sO], where the dipole-nucleon 
cross section is parametrized, in the Kharzeev-Levin (KL) approach it is the unintegrated gluon distribution which is 
parametrized. More details can be found in Ref . [25j . 

Another useful parametrization, which describes the HERA data, and therefore is valid for 10^"* < a; < 10^^ was 
constructed by Ivanov and Nikolaev (IN) [3^. We refer the reader for details to the original paper. 

In some of approaches one imposes the following relation between the standard coUinear distributions and UGDFs: 



(19) 



Due to its simplicity the Gaussian smearing of initial transverse momenta is a good and popular reference for other 
approaches. It allows to study phenomenologically the role of transverse momenta in several high-energy processes. 
We define a simple unintegrated gluon distribution: 



r. 



Gauss 



(x, k. 



xg^°^\x,nl) ■ fcauUkt) , 



(20) 



where g (x, Hp) is a standard coUinear (integrated) gluon distribution and fcaussiki) is a Gaussian two-dimensional 
function: 



f Gauss {kt ) 



1 



exp {-k'j2al) /tt . 



The UGDF defined by Eq.dSO]) and dUl) is normahzcd such that 



At small values of x the unintegrated gluon distribution can be obtained from integrated distribution as pi 

d{xg{x,fi'^)) 



T{x,kt) 



dn^ 



(21) 



(22) 



(23) 



This method cannot be directly used at small transverse momenta (small factorization scales) and must be supple- 
mented by a further prescription. One possible prescription is freezing of the gluon distribution at /c| < /i^^, another 
is a shift of the scale: fi^ n'^ + Of course /i^^ and /i^ are bigger than the lowest possible scale for standard 
coUinear distributions. This method cannot be also applied at larger x as here the scalling violation reverses and 
negative values are obtained. 

At intermediate and large x more careful methods must be used. Kwieciiiski has shown that the evolution equations 
for unintegrated parton distributions take a particularly simple form in the variable conjugated to the parton transverse 
momentum. In the impact-parameter space, the Kwiecihski equation takes the following simple form [ssj 









27rii^ Jo 




- fNs{x,b,fi^ 


dfs{x,b,fi'^) 


"s(m2) 




27r^2 








«s(/i^) f' 




27r/i2 y„ 



dzPqq{z) 



e{z - x) Jo((l - z)^ib) Jns ( -,b, 

V z 



dz{Q{z - x) Jo({l ~ z)iib) 



Pn 



[zP,,(z)-|-zPg,(z)]/5(a;,6,M') 



(24) 



dz< 9(z — x) Jo((l — z)/i6) 



P, 



gg 



+ P,,(^)/g(^,6,m') 
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We have introduced here the short-hand notation 

/nS = fu — fu, fd — fd J 

(25) 

fs = fu + fu + fd + fd + fs + fs ■ 

The unintegrated parton distributions in the impact factor representation are related to the famihar colhnear distri- 
butions as foUows 

Mx,b^O,fi^)^^Pk{x,^i^) . (26) 
On the other hand, the transverse momentum dependent UPDFs are related to the integrated parton distributions as 

a;pfc(x,/i^) = / dkf fk{x,kf,fX^) . (27) 
Jo 

The two possible representations are interrelated via Fourier-Bessel transform 

/"OO 

fk{x,k1,i?) = I dbbJo{ktb)fk{x,b,fi^) , 

I'oo (28) 
fk{x,b,i?)^ I dktktJo{ktb)fk{x,k^,fi^). 



Jo 

The index k above numerates either gluons (fc=0), quarks (fc > 0) or antiquarks (fc < 0). 

The perturbative solutions f^.^^^^ {x,b, do not include nonperturbative effects such as, for instance, intrinsic 
transverse momenta of partons in colliding hadrons. One of the reasons is e.g. internal motion of constituents of the 
proton. In order to include such effects we modify the perturbative solution /|"^''*(a;, 6, /i|,) and write the modified 

parton distributions fg{x,b,^'jp) in the simple factorized form 

fg{x,b,f,l) = fP^^\x,b,f^%) ■ F^P{b) . (29) 
In the present study we shall use the following functional form for the nonperturbative form factor 

^^rW=^""W-exp(^-^) • (30) 

In Eg. ([50)1 bo is the only free parameter. 

While physically fk{x, /i^) should be positive, there is no obvious reason for such a limitation for fk{x, b, fj?). 

In the following we use leading-order parton distributions from Ref. 34] as the initial condition for QCD evolution. 
The set of integro-differential equations in b-space was solved by the method based on the discretisation made with 
the help of the Chebyshev polynomials (see [33|). Then the unintegrated parton distributions were put on a grid in 
X, b and /z^ and the grid was used in practical applications for Chebyshev interpolation. 

For the calculation of the direct J/tp production here the parton distributions in momentum space are more useful. 
This calculation requires a time-consuming multi-dimensional integration. An explicit calculation of the Kwiecihski 
UPDFs via Fourier transform for needed in the main calculation values of {xi,kl f.) and (x2, A;| f) (see next section) 
is not possible. Therefore auxiliary grids of the momentum-representation UPDFs are prepared before the actual 
calculation of the cross sections. These grids are then used via a two-dimensional interpolation in the spaces (a:i, fcj () 
and (x2, fc| j) associated with each of the two incoming partons. 

The Kwieciriski unintegrated parton distributions were used recently in applications to cc photoproduction [SS*! , cc 
correlations in nucleon-nucleon collisions 25j , production of gauge bosons 36], production of Standard Model Higgs 
boson [33|, inclusive production of pions [3l|, production of direct photons [39]. Good agreement with experimental 
data was obtained in the case when the data exist. 

In the approach of Ref. [1^, based on leading-order perturbative solution of BFKL equation, the unintegrated gluon 
density J-g{x, k^ , ff^) is calculated as a convolution of the ordinary (colhnear) gluon density g{x,fi^) with universal 
weight factors: 

.f,(x,fc2,M')- f'g{v,klfi')-g{-,fi')drj, (31) 
Jx V V 

Gin, kl m') = ^ Jo(2w'a,ln(l/,7)ln(/x2/A:?)), fc? < n\ (32) 
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g{v,klf,^) = -^/o(2y'a,ln(l/,7)ln(fc2/^2))^ ^2 > ^2^ (33) 

where Jq and Iq stand for Bessel functions (of real and imaginary arguments, respectively), and as = ag/Sn. The LO 
GRV set 40] was used in our calculations as the input coUinear density. Here the value of as and the scale fj? are 
parameters of the model. The resulting unintegrated gluon distributions depend on them rather strongly. Sometimes 
for brevity we shall denote the distribution from Ref.[29|] by JB. 

III. RESULTS 

Now we shall compare contributions of different processes discussed in the previous section. Here a Monte Carlo 
method based on the VEGAS routine [4l| is used to allow an easy comparison of processes with different number of 
particles in the final state. In Fig. [3] we show contribution of different mechanisms discussed above to the rapidity 
distributions of J/tp meson for the RHIC energy. This calculation is based on so-called derivative UGDFs, i.e. the 
ones obtained by differentiating the standard coUinear distributions (see the previous section). 

In Fig[3]we show corresponding contributions to the transverse momentum distribution of the J/ip meson. In this 
exploratory calculation the cross section is integrated over the full range of rapidity. We obtain a rather surprising 
result that the sequential production of J/ip mesons via radiative decays of Xc mesons is comparable to or even 
dominates over the direct color-singlet contribution almost in the whole phase space. The reason can be seen in the 
fact that the production of Xcj states refers to much lower values of the final state invariant mass, << (p^ +pg)'^, 
giving emphasis to the small x region, where the gluon distributions are growing up. This property becomes even 
more pronounced as the 'direct' matrix element ([T]) vanishes when the emitted final gluon is soft. Our conclusion on 
the relative size of the direct and indirect contributions is compatible with the preliminary estimates obtained by the 
CDF collaboration 

We wish to note now some difficulties of the standard coUinear approach for the Xc mesons. The leading order 
contribution coming from the subprocess ([3|) shows unphysical ^-like pT spectrum. The usual excuse that the particles 
produced at zero pt disappear in the beam pipe and remain invisible does not work, because the decay products do 
have nonzero px and can be detected. At the same time, introducing the next-to-leading contributions (i.e., the 
processes with extra gluons in the final state) causes a problem of infrared divergences, which need artificial tricks to 
regularize them. 

It is well known that a large fraction of the -0' mesons decays into channels with J/ip (BR — 0.56 [13] )■ This 
contribution was not considered in the literature and requires a separate discussion. The inclusive cross section for ip' 
can be calculated in exactly the same way (color-singlet model) as the cross section for direct J/ tfj meson production. 
The decays of i/j' J/4' + X change the kinematics only slightly. Finally the ip' contribution constitutes about 25 % 
of the direct (color-singlet) production. 

Also the B-meson decay mechanism gives a sizeable contribution at large transverse momenta. 

Summarizing, at the RHIC energy the dominant production mechanisms are radiative decays of Xc(2^) and direct 
color-singlet mechanism. In the following we shall concentrate exclusively on these two dominant mechanisms. 

Let us start with color-singlet mechanism. In Fig|5] we present distributions in rapidity of J/ip produced by direct 
color-singlet mechanism for different UGDFs. The distribution obtained with Ivanov-Nikolaev (IN) glue exceeds the 
experimental PHENIX data [T§] , while the other theoretical distributions are smaller than experimental data. This is 
rather natural as contributions of other mechanisms are not included. The corresponding distributions in transverse 
momentum are shown in Fig[6]for two different intervals in rapidity. Very similar distributions are obtained for mid- 
and intermediate rapidity intervals. The result with Ivanov-Nikolaev UGDF exceeds the experimental data in the 
region of small transverse momenta. This is probably due to an extra nonperturbative contribution at small gluon 
transverse momenta (s^ . 

Now we shall show results obtained with different UGDFs for radiative decays of Xc(2^). The rapidity distribution 
of corresponding J/ip is shown in Figl?! Different UGDFs give a similar result. The distributions obtained with 
Ivanov-Nikolaev UGDF is slightly higher than those obtained with other distributions. In Figl5]we show distributions 
in transverse momentum of radiatively produced J/ip. The differences in the results for different UGDFs are up to a 
factor 2 or even larger. Again Ivanov-Nikolaev UGDF gives the highest cross section for small transverse momenta. 
The Bluemlein UGDF shown intentionally for large value of as = 0.6 (solid grey, green on line) gives completely 
wrong shape. The shape in this case depends strongly on the value of It would be much better for smaller values 
of as. 

Finally we would like to show how the sum of the two dominant contributions (direct color-singlet and radiative 
Xc(2^) decay) compares with the experimental data from RHIC. The distribution in rapidity is shown in Figl^land 
distributions in transverse momentum in Fig llOl The theoretical cross sections obtained with the Kwiecihski, BFKL 
and Kharzeev-Levin UGDFs stay slightly below the experimental data. This seems to be consistent with the fact 
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that the smaller contributions discussed in Fig[3] and Fig|4] are not included here. They are expected to produce 
contributions of the order of 20-30 % (see Fie:s l3l4p . 

At the RHIC energy W = 200 GeV the longitudinal momentum fractions of the order x ^ 10"^ - 10^^ come 
into game. This is the place where application of many UGDFs may be questionable. Let us concentrate now on 
Kwiecihski parton distributions, which are constructed for the region of x under discussion. In the left panel of 
Fig llll we show invariant cross section for the direct component as a function of J/ip transverse momentum pt for 
mid rapidity range -0.35 < y < 0.35. We show results for different factorization scales: fj? = 10 GeV^ (solid) and 
= 100 GeV^ (dashed). In the right panel of Fig[Tl]we show similar result for J/tp coming from the decays of the 
Xc(2^). Here the result depends more strongly on the choice of the scale. The solid line here corresponds to running 
factorization scale: /z^ = = ^^2+) +Pt- 

In Fig[T2] we compare the sum of both processes calculated with running factorization scale with the PHENIX 
experimental data. The calculation underestimates the data at small transverse momenta. This is most probably due 
to the omission of other components, especially the ■(/''-decay component. 

Let us concentrate now on the region of large transverse momenta of J/tp- In Fig |13l we show the contribution 
of J/4> from decays of the B and B mesons. The cross section for the bb is obtained with the Kwiecihski UGDF 
(fixed factorization and renormalization scales, /i^ = 4m^) within the fc^-factorization approach. The details of the 
calculation can be found in Ref.jl^. In the present illustrative calculation we neglect hadronization, i.e. we assume 
that the distribution of B (B) mesons is the same as the distribution of b (b) quarks. This seems justified for heavy 
quark to heavy meson transitions. There are several decay channels with final state J/ip. The inclusive branching 
ratio is known experimentally BR = 1.09% (20| . However, the momentum distribution of J/tp in the B meson center- 
of-mass system was not yet measured [sij . Here, in order to demonstrate the dependence on the details of the decay, 
we consider 3 academic models of the decays: (a) uniform distribution in p* (momentum of J/ip in the meson rest 
frame) in the interval {0,pmax) - dashed line, (b) uniform distribution inside the sphere with r&dius Pmax - dotted line, 
(c) distribution on the sphere with radius Pmax - dash-dotted line. Here Pmax is the momentum obtained assuming 
a two-body decay: B J/tpX. We assume the efi'ective mass of the state X to be mx = 0.5 GeV. As can be seen 
from the figure the B decays become an important ingredient at larger transverse momenta. There is relatively mild 
dependence on the details of the decay. However, these details may become important with a better statistics, when 
J/ip with Pt > 10 GeV will be measured. The present estimate of the B-decay contribution may be an underestimation 
because of the two following reasons: (a) it is based on leading-order approach, (b) choice of the renormalization scale 
(see above) [52j. Therefore at presently measured maximal transverse momenta oi J/tp pt ^ 8 GeV the B-decay 
contribution at the level of 20-30 % is not excluded. 

Let us concentrate now on correlations between produced J/tp and associated gluon(s). In Fig |14l we present two- 
dimensional distribution in transverse momentum of J/tp (pu) and transverse momentum of the associated (the gluon 
related with the matrix element) gluon {p2t) for two different scales of the Kwiecihski UGDF. The bigger the scale is, 
the bigger is the spread of the cross section in the [pit,P2t) space. This can be understood by the fact that the bigger 
scales means more gluonic emissions which statistically means the bigger spread. This figure is rather of academic 
value as in practice there are also gluons emitted in the process of the ladder-type emissions. Strictly speaking, the 
latter have to be described using a full gluon evolution generator. On the other hand, the relevant effects can also 
be estimated in an approximate way, as follows. On the average, the gluon transverse momentum increases from the 
proton line towards the hard interaction block (although there is no strict ordering in the transverse momentum in 
BFKL equation). So, it is most likely, that the last gluon in the parton ladder has the largest kt value. As a rough 
approximation, one can neglect the transverse momenta of all the other emitted gluons (note that the evolution is in 
the log(fct) space rather than kt space) and use the conservation law in the last splitting vertex to set the k't of the 
emitted gluon opposite to the kt of the gluon entering the partonic matrix element: k't ~ —kt- The latter is known 
from the unintegrated gluon distribution. This trick gives an estimate for the transverse momentum of the final state 
gluon jet. 

In Fig |15l we show distributions of the cross section on the plane pt{J /tp) x pt (matrix element gluon or last gluon in 
the ladder) for the Kwiecihski UGDF with running scale (left part) and BFKL UGDF (right part). Comparing these 
distributions we conclude that the gluons from the ladder (LFL - last from the ladder) contribute to lower transverse 
momenta than those associated with the matrix element g + g ^ J/tp + g (ME) for the Kwiecihski UGDF, where at 
Pf (gluon) > 5 GeV the matrix element gluons dominate over the ladder gluons. For the BFKL gluons the situation 
is much more complicated as here the distribution for ME gluons and LFL gluons are similar. 

In Fig[T6]we show average value of transverse momentum of the matrix element gluon (dashed line) and of the last 
gluon from the ladder (solid line) as a function of J/tp transverse momentum. These average values have completely 
different dependence on pt{J/tp). While average value of the LFL gluon transverse momentum is only weakly depen- 
dent on pt{ J/tp), the average value of the ME gluon transverse momentum grows monotonically with pt{J/tp). At low 
J/tp transverse momenta {pt{LFL)) ~ {pt{ME)). At higher J/tp transverse momenta {pt{LFL)) < {pt{ME)). For 
the Kwiecihski distribution this happens at smaller transverse momentum than for the BFKL UGDF. 
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Our calculations presented up to now show that the production of J /tp through radiative decays of Xc niesons is 
one of two dominant mechanisms. It would be worth to verify this theoretical prediction experimentally. This would 
require measuring the Xc mesons independently. In Fig ll7l we show the distributions in rapidity for Xc(0~''), Xc(l^) 
and Xc(2~''). These results were obtained with the Kwiecihski UGDF, which seems to be the most reliable for the 
RHIC energy range. 

For completeness in Fig llSI we show the corresponding distributions in transverse momentum. In this calculation 
-1 < y < 1. We wish to point out that the cross sections show no singularity at small transverse momentum. 
This contrasts with the collinear factorization predictions, which are either unphysicsl ((5-like) or divergent (if based 
on a 2 — > 2 subprocess g + g Xc + g)- There is also a significant difference in shape between the transverse 
momentum distribution for Xc(l^) meson and those for Xc(O^) and Xc(2~'") mesons. This property emerges from 
Landay-Yang theorem which prohibits the coupling of vector states to massless photons (just because of quantum 
numbers incompatible with Bose statistics). The production of Xc(l"'') states at small px is strongly suppressed 
because the initial gluons are almost on shell. The suppression goes away at higher pT, as the off-shellness of the 
initial gluons becomes larger. These features are discussed in detail in Ref. (43l |. 

In contrast to transverse momentum distribution of J/tp from the color-singlet mechanism, the distributions of Xc 
mesons (and consequently the distribution of Jf-ip from radiative decays) strongly depend on the model of UGDF. In 
particular, in the limiting case of vanishing initial gluon transverse momenta: da/cPpt oc 5'^[pt). For illustrating the 
effect quantitatively in FiglTHlwe present transverse momentum distributions of Xc(2^) for the Gaussian UGDF with 
different values of the smearing parameter ctq = 0.5, 1, 2 GeV. The example clearly demonstrates that a measurement 
of transverse momentum distribution of Xc mesons would open a new and unique possibility to test model unintegrated 
gluon distributions. 

In principle the Xc mesons (mainly Xc(l^) and Xc(2''')) can be identified via photon- J/^ decay channel. At RHIC 
the Xc production mechanism could be also identified using the tt+tt" and K^K^ final channels. The corresponding 
branching ratios are poj : 

Bi?(xc(0+) -> TT+TT") = 7.2 ± 0.6 X 10-3, BR{xc{0+) K+R-) = 5.4 ± 0.6 x lO'^, 
BR{xc{'2+) vr+vr-) = 2.14 ± 0.25 x lO^^, BR{xc{2+) K+R-) = 7.7 ± 1.4 x 10"''. 

Now we are coming to the issue of J/ip spin alignment, which was, and still is under intense debates in the literature. 
We want to stress once again that measuring the polarizaton of quarkonium states produced at high energies may 
serve as a crucial test discriminating the different concepts of parton dynamics. 

The polarization state of a vector meson is characterized by the spin alignment parameter a which is defined as a 
function of any kinematic variable as 

a{V) = [da/dV - Mgl / dV) / {da / dV + dcjL/dV), (34) 

where a is the reaction cross section, V \s a. selected kinematical variable and ctl is the part of cross section corre- 
sponding to mesons with longitudinal polarization (zero helicity state). The limiting values a = 1 and a = — 1 refer 
to the totally transverse and totally longitudinal polarizations. Here we consider only the behavior of a as a function 
of the J/'0 transverse momentum: V = |pt|- The experimental definition of a is based on measuring the angular 
distributions of the decay leptons 

dr{J/ij-^H+p,- )/dcos9 ^l + a cos^ 6, (35) 

where is the polar angle of the final state muon measured in the decaying meson rest frame. 

The results of our calculations for the kinematic conditions of RHIC are displayed in Figl^Ol In order to show the 
theoretical uncertainty band connected with the choice of UGDF, we use two different parametrizations, which are 
known to show the largest difference with each other, namely, the ones proposed in Refs. [l^ (called 'derivative' for 
brevity) and the one from Ref.[29j. 

The upper panel in Fig l20l shows the behavior of the spin alignment parameter a for J/ip mesons produced in 
the direct subprocess ([l} . The increase in the fraction of longitudinally polarised mesons comes from the increasing 
virtuality (and longitudinal polarization) of the initial gluons. These predictions shown here are also valid for ip' 
mesons. 

As far as the contribution from P-waves is concerned, nothing is known on the polarisation properties of their 
decays. If we assume that the quark spin is conserved in radiative transitions, and the emission of a photon only 
changes the quark orbital momentum (as it is known to be true in the electric dipole transitions in atomic physics, 
= 0, AL = ±1), then the predictions on a appear to be similar to those made for the direct channel (see 
lower panel in Figl^Ol dotted curves). If, on the contrary, we assume that the transition Xc^J/4'+J leads to a 
complete depolarization, then we arrive at a more moderate behavior of the parameter a (dash-dotted curves in 
Fig|20)l. The overall polarization remains slightly longitudinal (a ~ —0.1) in the whole range of due to the 'direct' 
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contribution. A comparison between the data on J/ij) and ■0' polarization at the Tevatron [8| seems to give support 
to the depolarization hypothesis. The difference between the J /ip and V'' polarization data can be naturally explained 
by the presence of the depolarizing contribution in the case of J /ip and the absence of this contribution in the case of 




FIG. 3: Contributions of different mechanisms for the production of J /ijj in da/dy distributions. In this calculation we have 
used simple "derivative UGDF" . 



IV. DISCUSSION AND CONCLUSIONS 



We have considered different mechanisms contributing to the inclusive production of J/tJj mesons in pp collisions 
at RHIC kinematics. The outcome of our study is the following. 

We have inspected the hierarchy of contributions and found that the dominant contribution to the cross section 
comes from radiative decay of Xc mesons, mainly from Xc(2^^) state. The second most important mechanism is the 
direct color-singlet production. The sequential process through the intermediate ip' turned out to be nonnegligible 
and constitutes about a quarter of the direct color-singlet contribution. To our knowledge, these processes were not 
included in previous calculations in the literature on the subject. 

As a by-product, we have demonstrated the advantage of the fct-factorization approach in calculating the Xc spectra: 
the latter can hardly be calculated in a cocsistent way in the collinear scheme. In order to verify the production 
mechanism suggested in our analysis, we have proposed an independent measurement of inclusive Xc cross sections in 
the TT+TT^ and K^K^ decay channels. 



12 




2 4 6 8 

~Pt (GeV) 



FIG. 4: Contributions of different mechanisms for the production of J /xp in da/dpt distributions. In this calculation we have 
used simple "derivative UGDF" . The cross section is integrated over the full range of rapidity. 

We have applied our approach to describe the data on inclusive J/V' production recently collected by the PHENIX 
Collaboration at the BNL. Both rapidity and transverse momentum distributions have been discussed. The new 
precise data at small J/ij} transverse momenta appeared to show very strong analysing power, imposing stringent 
constraints on unintegrated gluon distributions. The best description of the data is obtained with the UGDF proposed 
by Kwiecinski. 

Another piece of important information on the underlying gluon dynamics can be extracted from studying kinematic 
correlations between J/i/i mesons and coproduced gluon jets. In this paper we have presented our predictions for the 
two dominant contributing mechanisms. 

Finally, we have presented our predictions on the J/ij) spin alignment. The latter can serve as important test 
discriminating two different concepts of parton model. 

In the present paper we have discussed mechanisms of J /ip production in elementary collisions. We believe that 
our findings here may be also useful for nuclear collisions, where J/V' suppression was originally suggested as a useful 
indication of the presence of the quark-gluon plasma. 
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experimental data on inclusive J /ip production measured at RHIC. Thanks go to Jerzy Bartke for pointing to us the 
possibility that the ip'-decay contribution may be sizeable. This work was partially supported by the grant of the 
Polish Ministry of Scientific Research and Information Technology number 1 P03B 028 28. 
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FIG. 5: Direct color-singlet contribution to rapidity distribution of J/ip for different models of UGDFs. The solid (red on 
line) curve corresponds to the Kwiecihski UGDF, the dashed line to the Kharzeev-Levin UGDF, the dotted line to the BFKL 
UGDF, the dash-dotted line to the Ivanov-Nikolaev UGDF and the grey solid (green on line) curve to the Bluemlein UGDF. 
The V'' contribution is not included here. The new PHENIX data are shown as full circles. 
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the same as in Fig. [5] The tli' contribution is not included here. The new PHENIX data are shown as full circles. 
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FIG. 8: Xc-decay contribution to transverse momentum distribution of J /tj} for different models of UGDFs for different intervals 
in rapidity: (a) -0.35 < y < 0.35 (left panel), (b) 1.2 < \y\ < 2.2 (right panel). The meaning of the curves is the same as in 
Fig. [5] The new PHENIX data are shown as full circles. 
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FIG. 9: Direct color-singlet and Xc(2'^) contributions to rapidity distribution of J/ip for different models of UGDFs. The 
meaning of the curves is the same as in Fig. (5] The new PHENIX data are shown as full circles. 
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FIG. 10: Direct and Xc-decay contributions to transverse momentum distribution of J/ijj for different models of UGDFs for 
different intervals in rapidity: (a) -0.35 < y < 0.35 (left panel), (b) 1.2 < \y\ < 2.2 (right panel). The meaning of the curves is 
the same as in Fig. [5] The new PHENIX data are shown as full circles. 
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FIG. 11: Factorization scale dependence of the transverse momentum distribution for Kwiecinski UGDF. The mid rapidity 
range -0.35 < y < 0.35 was taken as an example. The left panel is for direct production and the right panel for the Xc{'2'^) 
decay mechanism. The solid and dashed curves are for fi^ = 10 GeV'^ and for fi^ = 100 GeV^, respectively. In this calculation 
bo = 1 GeV-\ 
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FIG. 12: Invariant cross section for the Kwieciriski UGDF with running scale. The left panel is for central rapidity range (-0.35 
< y < 0.35) and the right panel for intermediate rapidity range (1.2 < y < 2.2). The direct contribution is denoted by the 
dashed line, the Xc(2''')-decay contribution by the dotted line and the sum of both by the solid line. 
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FIG. 13: Invariant cross section for J/?/> from decays of the B mesons as a function of pt for midrapidity and intermediate 
rapidity intervals. Kwiecinski UGDFs are used with factorization scale fj,^ = ^ml. Different decay models are described in the 
text. 
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FIG. 15: Two-dimensional distribution of the J/^ and gluon transverse momenta. The left-top panel is for Kwieciriski UGDF 
(running scale) and matrix element gluon, the right-top panel is for the BFKL UGDF and matrix element gluon, the left-bottom 
panel is for the Kwiecinski UGDF (running scale) and "last from the ladder" gluon and the right-bottom panel is for the BFKL 
UGDF and "last from the ladder" gluon. 
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FIG. 16: Average transverse momenta of the ME (dashed) and LFL (sohd) gluons as a function of the J/t/j transverse momentum 
for the Kwiecmski UGDF with running scale (left panel) and the BFKL UGDF (right panel) at the RHIC energy W = 200 
GeV. 
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FIG. 17: Rapidity distribution of Xc(0+) (dashed), Xc(l+) (dotted) and Xc(2+) (solid) for the RHIC energy obtained with the 
Kwiecinski UGDF (60 = 1 GeV-\ fJ.^ = PtiXc))- 
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FIG. 18: Transverse momentum distribution of Xc(0"'') (dashed), Xci^^) (dotted) and Xc(2"'") (solid) for the RHIC energy 
obtained with the Kwiecihski UGDF (6o = 1 GeV"\ =P?(Xc))- 
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FIG. 19: Transverse momentum distribution of Xc{2'^) for the RHIC energy obtained with the Gaussian UGDF and different 
values of ao = 0.5, 1.0, 2.0 GeV. 
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FIG. 20: Predictions for the spin-alignment parameter a for J/-)/; and W = 200 GeV. Thick lines correspond to the Bluemlein 
parametrization [2^ and the thin lines correspond to the derivative UGDF parametrization and the GRV coUinear distribution. 
The top panel is for direct contribution only. The bottom panel includes the feed-down from Xc decays taken into account. The 
dotted lines are for the quark spin conservation hypothesis, and the dash-dotted lines are for the full depolarization hypothesis. 



